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Abstract 

The paper proposes a control-theoretic framework for verification of numerical software systems, and puts 
forward software verification as an important application of control and systems theory. The idea is to transfer 
Lyapunov functions and the associated computational techniques from control systems analysis and convex op- 
timization to verification of various software safety and performance specifications. These include but are not 
limited to absence of overflow, absence of division-by-zero, termination in finite time, presence of dead-code, and 
certain user-specified assertions. Central to this framework are Lyapunov invariants. These are properly constructed 
functions of the program variables, and satisfy certain properties — resembling those of Lyapunov functions — along 
the execution trace. The search for the invariants can be formulated as a convex optimization problem. If the 
associated optimization problem is feasible, the result is a certificate for the specification. 

Index Terms 

Software Verification, Lyapunov Invariants, Convex Optimization. 

I. Introduction 

SOFTWARE in safety-critical systems implement complex algorithms and feedback laws that control 
the interaction of physical devices with their environments. Examples of such systems are abundant 
in aerospace, automotive, and medical applications. The range of theoretical and practical issues that arise 
in analysis, design, and implementation of safety-critical software systems is extensive, see, e.g., OTA . 
P4ll , and Il27ll . While safety-critical software must satisfy various resource allocation, timing, scheduling, 
and fault tolerance constraints, the foremost requirement is that it must be free of run-time errors. 

A. Overview of Existing Methods 

1) Formal Methods: Formal verification methods are model-based techniques [fSTTl . Il48l . iflTTl for 

proving or disproving that a mathematical model of a software (or hardware) satisfies a given specification, 
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i.e., a mathematical expression of a desired behavior. The approach adopted in this paper too, falls under 
this category. Herein, we briefly review model checking and abstract interpretation. 

a) Model Checking: In model checking [15] the system is modeled as a finite state transition system 
and the specifications are expressed in some form of logic formulae, e.g., temporal or propositional logic. 
The verification problem then reduces to a graph search, and symbolic algorithms are used to perform 
an exhaustive exploration of all possible states. Model checking has proven to be a powerful technique 
for verification of circuits [fT4|| - security and communication protocols ll38l . Il45ll and stochastic processes 
flU. Nevertheless, when the program has non-integer variables, or when the state space is continuous, 
model checking is not directly applicable. In such cases, combinations of various abstraction techniques 
and model checking have been proposed (31, [fT9ll , I162R : scalability, however, remains a challenge. 

b) Abstract Interpretation: is a theory for formal approximation of the operational semantics of 
computer programs in a systematic way ifToll . Construction of abstract models involves abstraction of 
domains — typically in the form of a combination of sign, interval, polyhedral, and congruence abstractions 
of sets of data — and functions. A system of fixed-point equations is then generated by symbolic for- 
ward/backward executions of the abstract model. An iterative equation solving procedure, e.g., Newton's 
method, is used for solving the nonlinear system of equations, the solution of which results in an inductive 
invariant assertion, which is then used for checking the specifications. In practice, to guarantee finite 
convergence of the iterates, narrowing (outer approximation) operators are constructed to estimate the 
solution, followed by widening (inner approximation) to improve the estimate IfTTll . This compromise can 
be a source of conservatism in analysis Q. Nevertheless, these methods have been used in practice for 
verification of limited properties of embedded software of commercial aircraft ||8l. 

Alternative formal methods can be found in the computer science literature mostly under deductive 
verification 071 . type inference 11521 . and dataflow analysis 11281 . These methods share extensive similar- 
ities in that a notion of program abstraction and symbolic execution or constraint propagation is present 
in all of them. Further details and discussions of the methodologies can be found in ifTTl . and 11481 . 
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2) System Theoretic Methods: While software analysis has been the subject of an extensive body of 
research in computer science, treatment of the topic in the control systems community has been less 
systematic. The relevant results in the systems and control literature can be found in the field of hybrid 
systems [12J. Much of the available techniques for safety verification of hybrid systems are explicitly 
or implicitly based on computation of the reachable sets, either exactly or approximately. These include 
but are not limited to techniques based on quantifier elimination 041 . ellipsoidal calculus [|32l . and 
mathematical programming [6]. Alternative approaches aim at establishing properties of hybrid systems 
through barrier certificates [|53l , numerical computation of Lyapunov functions [fTTI . [|29l , or by combined 
use of bi simulation mechanisms and Lyapunov techniques ll24l . [|33l . [|62l . J3J. 

Inspired by the concept of Lyapunov functions in stability analysis of nonlinear dynamical systems 
ll30ll . in this paper we propose Lyapunov invariants for analysis of computer programs. While Lyapunov 
functions and similar concepts have been used in verification of stability or temporal properties of system 
level descriptions of hybrid systems ll55ll . ifTTTl . E9ll . to the best of our knowledge, this paper is the first 
to present a systematic framework based on Lyapunov invariance and convex optimization for verification 
of a broad range of code-level specifications for computer programs^] Accordingly, it is in the systematic 
integration of new ideas and some well-known tools within a unified software analysis framework that we 
see the main contribution of our work, and not in carrying through the proofs of the underlying theorems 
and propositions. The introduction and development of such framework provides an opportunity for the 
field of control to systematically address a problem of great practical significance and interest to both 
computer science and engineering communities. The framework can be summarized as follows: 

1) Dynamical system interpretation and modeling (Section [El]). We introduce generic dynamical system 
representations of programs, along with specific modeling languages which include Mixed-Integer 
Linear Models (MILM), Graph Models, and MIL-over- Graph Hybrid Models (MIL-GHM). 



2) Lyapunov invariants as behavior certificates for computer programs (Section III). Analogous to a 
Lyapunov function, a Lyapunov invariant is a real-valued function of the program variables, and 
satisfies a difference inequality along the trace of the program. It is shown that such functions can 

'This paper constitutes the synthesis and extension of ideas and computational techniques expressed in the workshop |20|, and subsequent 
conference papers |56|, |57) and 1581 . Some of the ideas presented in |20|, |56), and 1571 were independently reported in II 81 . 
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be formulated for verification of various specifications. 



3) A computational procedure for finding the Lyapunov invariants (Section IV). The procedure is 
standard and constitutes these steps: (i) Restricting the search space to a linear subspace. (ii) Using 
convex relaxation techniques to formulate the search problem as a convex optimization problem, 
e.g., a Linear Program (LP) Q, Semidefmite Program (SDP) OH, j63l, or a Sum-of-Squares (SOS) 
program [|49l . (iii) Using convex optimization software for numerical computation of the certificates. 

II. Dynamical System Interpretation and Modeling of Computer Programs 
We interpret computer programs as discrete-time dynamical systems and introduce generic models that 

formalize this interpretation. We then introduce MILMs, Graph Models, and MIL-GHMs as structured 

cases of the generic models. The specific modeling languages are used for computational purposes. 

A. Generic Models 

1 ) Concrete Representation of Computer Programs: We will consider generic models defined by a 
finite state space set X with selected subsets X C X of initial states, and X^ C X of terminal states, 
and by a set-valued state transition function / 2 X , such that f(x) C X^^Wx G X^. We denote 

such dynamical systems by S(X, f, Xq.X^). 

Definition 1: The dynamical system S(X, f, X , Xoo) is a C-representation of a computer program 
V, if the set of all sequences that can be generated by V is equal to the set of all sequences X = 
(x(0), x(l), . . . , . . . ) of elements from X, satisfying 

x(0)ex cx, x (t + 1) e / (x (*)) Vt e z+ (l) 

The uncertainty in a;(0) allows for dependence of the program on different initial conditions, and the 
uncertainty in / models dependence on parameters, as well as the ability to respond to real-time inputs. 

Example 1: Integer Division (adopted from 115111 ): The functionality of Program 1 is to compute the 
result of the integer division of dd (dividend) by dr (divisor). A C -representation of the program is 
displayed alongside. Note that if dd > 0, and dr < 0, then the program never exits the "while" loop and 
the value of q keeps increasing, eventually leading to either an overflow or an erroneous answer. The 
program terminates if dd and dr are positive. 
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int IntegerDivision ( int dd, int dr 
{hit q = {0}; int r = {dd}; 
while (r >= dr) 

{ q = q+i; 

r = r — dr; } 
return r; } 



Z = Zn [-32768, 32767] 






X = z 4 






\ r {(AA A *. ^ v.\ r~ "V 1 ^ fl 

ao = {(dd, dr, q, r) G A | q = (J, i = 


ad} 




X x = {(dd,dr,q,r) G X | r < dr} 






f (dd, dr, q + 1, r 


- dr), 


(dd, dr, q, r) G X\a' 00 


/:(dd,dr,q,r)^| {dd ^ q ^ 




(dd, dr, q, r) G 



Program 1: The Integer Division Program (left) and its Dynamical System Model (right) 
2) Abstract Representation of Computer Programs: In a C-representation, the elements of the state 
space X belong to a finite subset of the set of rational numbers that can be represented by a fixed 
number of bits in a specific arithmetic framework, e.g., fixed-point or floating-point arithmetic. When 
the elements of X are non-integers, due to the quantization effects, the set-valued map / often defines 
very complicated dependencies between the elements of X, even for simple programs involving only 
elementary arithmetic operations. An abstract model over-approximates the behavior set in the interest of 
tractability. The drawbacks are conservatism of the analysis and (potentially) undecidability. Nevertheless, 
abstractions in the form of formal over- approximations make it possible to formulate computationally 
tractable, sufficient conditions for a verification problem that would otherwise be intractable. 

Definition 2: Given a program V and its C-representation <S(X, /, X , Xoo), we say that <S(X, /, X , Xoo) 
is an ^-representation, i.e., an abstraction of V, if X C X, X C X , and f(x) C f[x) for all x E X, 
and the following condition holds: 

Iconic^. (2) 



Thus, every trajectory of the actual program is also a trajectory of the abstract model. The definition 
of Xqo is slightly more subtle. For proving Finite-Time Termination (FTT), we need to be able to infer 
that if all the trajectories of S eventually enter !„, then all trajectories of S will eventually enter I^. 
It is tempting to require that I M C however, this may not be possible as I ro is often a discrete set, 
while Xqo is dense in the domain of real numbers. The definition of Xqo as in ((2} resolves this issue. 

Construction of <S(X, /, X , X^) from <S(X, /, X , X^) involves abstraction of each of the elements 
X, /, X , Xqo in a way that is consistent with Definition[2] Abstraction of the state space X often involves 
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replacing the domain of floats or integers or a combination of these by the domain of real numbers. 
Abstraction of X or Xoo often involves a combination of domain abstractions and abstraction of functions 
that define these sets. Semialgebraic set- valued abstractions of some commonly used nonlinearities is 
presented in Appendix I. Also, abstractions of fixed-point and floating point arithmetic operations based 
on ideas from |42| and ll43l are discussed in Appendix I. A case study in application of these methods 



to analysis of a program with floating-point operations is presented in Section V-B 
B. Specific Models of Computer Programs 

Specific modeling languages are particularly useful for automating the proof process in a computational 
framework. Here, three specific modeling languages are proposed: Mixed-Integer Linear Models (MILM), 
Graph Models, and Mixed-Integer Linear over Graph Hybrid Models (MIL-GHM). 

1) Mixed-Integer Linear Model (MILM): Proposing MILMs for software modeling and analysis is 

motivated by the observation that by imposing linear equality constraints on boolean and continuous 

variables over a quasi-hypercube, one can obtain a relatively compact representation of arbitrary piecewise 

affine functions defined over compact polytopic subsets of Euclidean spaces (Proposition^. The earliest 

reference to the statement of universality of MILMs appears to be [|46ll , in which a constructive proof is 

given for the one-dimensional case. A constructive proof for the general case is given in Appendix II. 
Proposition 1: Universality of Mixed-Integer Linear Models. Let / : X h-» W 1 be a piecewise affine 

map with a closed graph, defined on a compact state space X C [—1,1]™, consisting of a finite union of 

compact poly topes. That is: 

/ (x) G 2AiX + 2Bi subject to x E X h i G Z (1, N) 

where, each Xi is a compact polytopic set. Then, / can be specified precisely, by imposing linear 
equality constraints on a finite number of binary and continuous variables ranging over compact intervals. 
Specifically, there exist matrices F and H, such that the following two sets are equal: 

G 1 = {(x,f(x)) \xeX} 

G 2 = {(x,y) \F[ X w v if =y, H[ x w v if = 0, (w,v) G [-1, lp x {-1, 1}""} 
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Mixed Logical Dynamical Systems (MLDS) with similar structure were considered in [|5l for analysis 
of a class of hybrid systems. The main contribution here is in the application of the model to software 
analysis. A MIL model of a computer program is defined via the following elements: 

1) The state space X c [—1, l] n . 

2) Letting n e — n + n w + n v + 1, the state transition function / : X h-» 2 X is defined by two matrices 
F, and H of dimensions n-by-n e and n#-by-n e respectively, according to: 

f(x) E w v if | H[ x w v if = 0, (w,v) E [-1,1]"™ x {-l,l} n "|. (3) 

3) The set of initial conditions is defined via either of the following: 

a) If X is finite with a small cardinality, then it can be conveniently specified by extension. We 



will see in Section IV that per each element of X , one constraint needs to be included in 
the set of constraints of the optimization problem associated with the verification task, 
b) If X is not finite, or \X \ is too large, an abstraction of X can be specified by a matrix 
H E R UH o xrie which defines a union of compact polytopes in the following way: 

X = {x E X | H [x w v 1 ] T = 0, (w,v) E [-1,1]"" x {-1,1}""}. (4) 
4) The set of terminal states X^ is defined by 

X 00 = {xEX \ H[ x w v 1 f ^ 0, E [-1, lp , Vv E {-1, lp}. (5) 



Therefore, S(X, f, X^X^) is well defined. A compact description of a MILM of a program is either 
of the form S (F, H, H ,n,n w ,n v ) , or of the form S (F, H, X ,n,n w ,n v ). The MILMs can represent a 
broad range of computer programs of interest in control applications, including but not limited to control 
programs of gain scheduled linear systems in embedded applications. In addition, generalization of the 
model to programs with piecewise affine dynamics subject to quadratic constraints is straightforward. 



Example 2: A MILM of an abstraction of the IntegerDivision program (Program 1: Section II-A), with 
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all the integer variables replaced with real variables, is given by S (F, H, H , 4, 3, 0) , where 

H = 
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Here, M is a scaling parameter used for bringing all the variables within the interval [— 1, 1] . 
2) Graph Model: Practical considerations such as universality and strong resemblance to the natural 

flow of computer code render graph models an attractive and convenient model for software analysis. 

Before we proceed, for convenience, we introduce the following notation: P r (i, x) denotes the projection 

operator defined as P r (i, x) = x, for all % G ZU {n} , and all x G M. n . 

A graph model is defined on a directed graph G (Af, £) with the following elements: 

1) A set of nodes Af = {0} U {1, . . . , m} U {x} . These can be thought of as line numbers or code 
locations. Nodes and x are starting and terminal nodes, respectively. The only possible transition 

from node n is the identity transition to node x . 

2) A set of edges £ = k) \ i G Af, j G O (i)} , where the outgoing set O (i) is the set of all 

nodes to which transition from node % is possible in one step. Definition of the incoming set X {%) 
is analogous. The third element in the triplet k) is the index for the A;th edge between % and 

j, and Aji = {k | (i,j,k) £ £} . 

3) A set of program variables xi e C R, I e Z (1, n) . Given Af and n, the state space of a graph 

model is X = Af x Q n . The state x — (i, x) of a graph model has therefore, two components: The 

discrete component % e Af, and the continuous component x £ il n C M. n . 

4) A set of transition labels T ^ assigned to every edge (i, j, k) e £, where maps x to the set T^x = 



{7j. (x,iu,u) | G S^}, where (w,v) e [-1, l]" 1 " x {-1, l} n " , and : M n+n - 



i — y 



is a polynomial function and 5^ is a semialgebraic set. If T Vl is a deterministic map, we drop S: 



and define = 7j (x). 

5) A set of passport labels 11^ assigned to all edges k) e where 11^ is a semialgebraic set. A 
state transition along edge k) is possible if and only if x G n^. 
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6) Semialgebraic invariant sets Xj C fl n , % g J\f are assigned to every node on the graph, such that 

P r (i,x) G Xi. Equivalently, a state x = (i,x) satisfying x G X\Xi is unreachable. 
Therefore, a graph model is a well-defined specific case of the generic model S(X, f, XojXqo), with 
X = M x Q n , X = {0} x X 9 , l 00 = {N}xX M and/:1^2 x defined as: 

f(x) = f(i,x) = {(j,l* ji x) \jeO(i), xen^nXi}. (6) 

Conceptually similar models, namely control flow graphs, have been reported in [51] (and the references 
therein) for software verification, and in fl2), [fT3ll for modeling and verification of hybrid systems. 
Remarks 

- The invariant set of node contains all the available information about the initial conditions of the 

program variables: P r (0,x) G X$. 

- Multiple edges between nodes enable modeling of logical "or" or "xor" type conditional transitions. 

This allows for modeling systems with nondeterministic discrete transitions. 

- The transition label may represent a simple update rule which depends on the real-time input. 

For instance, if T = Ax + Bw, and S = R n x [—1, 1] , then x H- {Ax + Bw \ w G [—1, 1]} . 
In other cases, T 4 may represent an abstraction of a nonlinearity. For instance, the assignment 



x \-> sin (x) can be abstracted by x i-> {T (x, w) \ (x, w) G S} (see Eqn. (63) in Appendix I). 



- Graph models with state-dependent or time- varying transitions labels arise, for instance, in modeling 

computer programs that involve arithmetic operations with array elements. For instance, consider 

the case where the transition functions are parametrized, and the parameters are drawn from a finite 

set, e.g., a multidimensional array. The dependence on the array elements can be random, resulting 

time-varying transition labels, or it can be in the form of a functional dependence on other program 

variables. Systematic treatment of such models is discussed in Appendix I. 

Before we proceed, we introduce the following notation: Given a semialgebraic set IT, and a polynomial 

function r : R n ^ ]R n , we denote by II (r) , the set: il(r) = {x | r (x) G 11} . 

a) Construction of Simple Invariant Sets: Simple invariant sets can be included in the model if they 

are readily available or easily computable. Even trivial invariants can simplify the analysis and improve 
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the chances of finding stronger invariants via convex relaxations, e.g., the 5-Procedure (cf. Section IV). 

- Simple invariant sets may be provided by the programmer. These can be trivial sets representing 
simple algebraic relations between variables, or they can be more complicated relationships that 
reflect the programmer's knowledge about the functionality and behavior of the program. 

- Invariant Propagation: Assuming that are deterministic and invertible, the set 

x t = U n*((i*)- 1 ) (7) 

is an invariant set for node i. Furthermore, if the invariant sets Xj are strict subsets of £l n for all 
j el (i) , then ^ can be improved. Specifically, the set 

*i= U nj-ter^n^r^)- 1 ) (8) 

is an invariant set for node i. Note that it is sufficient that the restriction of T£ to the lower 
dimensional spaces in the domains of 11^ and Xj be invertible. 

- Preserving Equality Constraints: Simple assignments of the form Tb : x\ M- / (y m ) result in invariant 
sets of the form Xj = {x \ x t — f (y m ) = 0} at node i, provided that does not simultaneously 
update y m . Formally, let be such that {T^x)i—xi is non-zero for at most one element I £ Z (1, n) , 
and that {T^x)i is independent of X[. Then, the following set is an invariant set at node % : 

*i= U {*l [IS-/] s = o} 

jeX(i), k^Aij 

3) Mixed-Integer Linear over Graph Hybrid Model (MIL-GHM): The MIL-GHMs are graph models 
in which the effects of several lines and/or functions of code are compactly represented via a MILM. As 
a result, the graphs in such models have edges (possibly self-edges) that are labeled with matrices F and 
H corresponding to a MILM as the transition and passport labels. Such models combine the flexibility 
provided by graph models and the compactness of MILMs. An example is presented in Section |VJ 

C. Specifications 

The specification that can be verified in our framework can generically be described as unreachability 
and finite-time termination. 
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Definition 3: A Program V = S(X, f, X , X^) is said to satisfy the unreachability property with 

respect to a subset X_ C X, if for every trajectory X = x(-) of ([I]), and every t G Z+, does 

not belong to X_. A program V = <S(X, j^X^X^ is said to terminate in finite time if every solution 

X = x (•) of ([I]) satisfies G for some t G Z + . 

Several critical specifications associated with runtime errors are special cases of unreachability. 

1) Overflow: Absence of overflow can be characterized as a special case of unreachability by defining: 

X_ = [x e X \ HaT 1 ^! > 1, a — diagjaj}} 

where a- t > is the overflow limit for variable i. We say that a y specifies the overflow limit. 

2) Out-of-Bounds Array Indexing: An out-of-bounds array indexing error occurs when a variable 

exceeding the length of an array, references an element of the array. Assuming that xi is the corresponding 
integer index and L is the array length, one must verify that x t does not exceed L at location i, where 
referencing occurs. This can be accomplished by defining X_ = {(i,x) G X \ \xi\ > L} over a graph 

model and proving that X_ is unreachable. This is also similar to "assertion checking" defined next. 

3) Program Assertions: An assertion is a mathematical expression whose validity at a specific location 

in the code must be verified. It usually indicates the programmer's expectation from the behavior of the 
program. We consider assertions that are in the form of semialgebraic set memberships. Using graph 
models, this is done as follows: 

at location i : assert x G A4 define X_ = {(i,x) G X \ x G X\Ai} , 
at location i : assert x A{ =^> define X_ = {(i, x) G X \ x G Ai] . 
In particular, safety assertions for division-by-zero or taking the square root (or logarithm) of positive 
variables are standard and must be automatically included in numerical programs (cf. Sec. III-A Table [i]). 



4) Program Invariants: A program invariant is a property that holds throughout the execution of the 
program. The property indicates that the variables reside in a semialgebraic subset Xj C X. Essentially, 
any method that is used for verifying unreachability of a subset X_ C X, can be applied for verifying 
invariance of Xj by defining X_ = X\Xj, and vice versa. 
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D. Implications of the Abstractions 

For mathematical correctness, we must show that if an ^4-representation of a program satisfies the 

unreachability and FTT specifications, then so does the C-representation, i.e., the actual program. This is 

established in the following proposition. 

Proposition 2: Let S(X, f, X , X^) be an ^.-representation of program V with C-representation 

S(X, f, Xq, Xoo). Let X_ C X and X_ C X be such that X_ C X— Assume that the unreachability 

property w.r.t. X_ has been verified for S. Then, V satisfies the unreachability property w.r.t. X_. 

Moreover, if the FTT property holds for S, then V terminates in finite time. 

Proof: See Appendix EE. ■ 

Since we are not concerned with undecidability issues, and in light of Proposition [2} we will not 

differentiate between abstract or concrete representations in the remainder of this paper. 

III. Lyapunov Invariants as Behavior Certificates 

Analogous to a Lyapunov function, a Lyapunov invariant is a real-valued function of the program 

variables satisfying a difference inequality along the execution trace. 

Definition 4: A (9, p) -Lyapunov invariant for S(X, f, X , X^) is a function V:Xi->M such that 

V (x+) - 9V (x) < -p Vx G X, x + G f (x) : x £ X M . (9) 

where (9, fx) E [0, oo) x [0,oo). Thus, a Lyapunov invariant satisfies the difference inequality ^ along 
the trajectories of S until they reach a terminal state X^. 

It follows from Definition [4] that a Lyapunov invariant is not necessarily nonnegative, or bounded from 
below, and in general it need not be monotonically decreasing. While the zero level set of V defines an 
invariant set in the sense that V (xk) < implies V (xk+i) < 0, for all I > 0, monotonicity depends on 
9 and the initial condition. For instance, if V (x ) < 0, Vx G X 0l then ^ implies that V (x) < along 
the trajectories of S, however, V (x) may not be monotonic if 9 < 1, though it will be monotonic for 
9 > 1. Furthermore, the level sets of a Lyapunov invariant need not be bounded closed curves. 
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Proposition [3] (to follow) formalizes the interpretation of Definition [4] for the specific modeling lan- 
guages. Natural Lyapunov invariants for graph models are functions of the form 

V{x)=V (z, x) =<Ti{x), ie N, (10) 

which assign a polynomial Lyapunov function to every node i E M on the graph G (M, S) . 

Proposition 3: Let S (F, H, X 0l n, n W) n v ) and properly labeled graph G (Af, S) be the MIL and graph 

models for a computer program V. The function V : [—1, 1]" n- ]R is a (9, ^-Lyapunov invariant for V 

if it satisfies: 

V(Fx e ) - 9V (x) < -fi, V (x, x e ) e [-1, l] n x 5, 

where 

S = {(>, w,v,l) | H[ x w v 1 f = 0, (w, v) e [-1, lp x {-1, lp}. 



The function V : J\fxW. n i— >• M, satisfying (10) is a (0, /z)-Lyapunov invariant for P if 

(Tj(x+) - Boi (x) < -fi, V k) e £, (x, x+) G (Xi n Iljj x T^x. (11) 
Note that a generalization of ([9]) allows for 6 1 and [i to depend on the state x, although simultaneous 
search for 9 (x) and V (x) leads to non-convex conditions, unless the dependence of 9 on x is fixed 
a-priori. We allow for dependence of 9 on the discrete component of the state in the following way: 

<Tj(x+) - 9%a t (x) < -tip, V k) e £, (x,x + ) e (X t n njj x T^x (12) 
A. Behavior Certificates 

1) Finite-Time Termination (FTT) Certificates: The following proposition is applicable to FTT analysis 
of both finite and infinite state models. 

Proposition 4: Finite-Time Termination. Consider a program V, and its dynamical system model 

S(X, f, XojXqc). If there exists a (9, /i)-Lyapunov invariant V : X i— > K, uniformly bounded on X\Xoo, 

satisfying (|9]) and the following conditions 

V (x) < — r/ < 0, Vi6 X (13) 

^ + (^-i)II^IL>° ( 14 ) 

max (//, ?y) > (15) 
where HVH^ = sup V (x) < oo, then P terminates in finite time, and an upper-bound on the number 
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of iterations is given by 

f tog(Ai + (0-l)||V||J-]og(Ai) 



log 9 



9 ^ 1, /x > 



T = < 



MMJ^jgg (V) n,, „ n (16) 



Proof: The proof is presented in Appendix EL ■ 
When the state- space X is finite, or when the Lyapunov invariant V is only a function of a subset of 
the variables that assume values in a finite set, e.g., integer counters, it follows from Proposition [4] that V 
being a (6>, /i)-Lyapunov invariant for any > 1 and \i > is sufficient for certifying FTT, and uniform 
boundedness of V need not be established a-priori. 

Example 3: Consider the IntegerDivision program presented in Example [TJ The function V : X >->■ R, 
defined according to V : (dd, dr, q, r) i— > r is a (1, dr) -Lyapunov invariant for IntegerDivision: at every 
step, V decreases by dr > 0. Since X is finite, the program IntegerDivision terminates in finite time. This, 
however, only proves absence of infinite loops. The program could terminate with an overflow. 

2) Separating Manifolds and Certificates of Boundedness: Let V be a Lyapunov invariant satisfying 
^ with 9 = 1. The level sets of V, defined by C r (V) = {i 6 I : V(x) < r}, are invariant with respect 
to ([T]) in the sense that x(t + 1) G C r (V) whenever x(t) E C r (V). However, for r = 0, the level sets 
C r (V) remain invariant with respect to ([T]) for any nonnegative 9. This is an important property with 
the implication that 9 = 1 (i.e., monotonicity) is not necessary for establishing a separating manifold 
between the reachable set and the unsafe regions of the state space (cf. Theorem [T]). 

Theorem 1: Lyapunov Invariants as Separating Manifolds. Let V denote the set of all (9,p)- 
Lyapunov invariants satisfying (j9| for program V = S(X, /, X , Xco). Let / be the identity map, and for 
he {f,I} define 

h- 1 (X_) = {x G X | h (x) n X_ ^ 0} . 
A subset X_ C X, where X_ fl X = can never be reached along the trajectories of V, if there exists 
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V E V satisfying 



supV(x) < inf V(x) (17) 

x£X x&h^{X^) 



and either 9 = 1, or one of the following two conditions hold: 



(I) 9 < 1 and inf V(x) > 0. (18) 

xeh-^x.) 

(II) 9 > 1 and supV(x) < 0. (19) 

x£X 



Proof: The proof is presented in Appendix EL ■ 
The following corollary is based on Theorem [T] and Proposition [4] and presents computationally 



implementable criteria (cf. Section IV) for simultaneously establishing FTT and absence of overflow. 

Corollary 1: Overflow and FTT Analysis Consider a program V, and its dynamical system model 
S(X, f, X , Xoo). Let a > be a diagonal matrix specifying the overflow limit, and let X_ = {x E 
x I ll a ~ 1;E lloo > !}• Let Q e NU{oo} , h E {/, /} , and let the function 7:li->lbea (0,/x)-Lyapunov 
invariant for 5 satisfying 

V (x) < Vx E X Q . (20) 

> supjUa" 1 ^^)^ - l} VxeX (21) 

Then, an overflow runtime error will not occur during any execution of V. In addition, if /i > and 

ji + 9 > 1, then, "P terminates in at most T u iterations where T u = fi~ l if 9 = 1, and for 6 1 ^ 1 we have: 

T = log (//+(g-l)||y||J-log// < log (// + g-l)-log/x 
log 6* ~~ log6> 

where || V || ^ = sup 

xex\{x_ux oc } 

Proof: The proof is presented in Appendix II. ■ 
Remarks 

- Application of Corollary [T] with h = f typically leads to much less conservative results compared 
with h = I, though the computational costs are also higher. 
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- An alternative, less conservative formulation is obtained when (21) is replaced by n constraints in 
the following way: 

V(x) > sup{|«~ 1 /i(x) fc | - 1} Vx G X, k G Z(l,n) (23) 
where a k is the overflow limit for scalar variable x k . Since a k x x k is a scalar quantity, \\a k 1 h (x) k \\ q = 



a k h(x) k \ and (23) can be used in conjunction with quadratic Lyapunov invariants and the 



5-Procedure for convex relaxation to formulate the problem of searching for V as semidefinite 



optimization problem (cf. Section IV). However, (23) is computationally more expensive than (21 ), 



as it increases the number of constraints. 
- An even less conservative but computationally more demanding alternative would be to construct 
a different Lyapunov invariant for each variable and rewrite (23), (21), (20) along with ([9]), for n 
different functions V k (.). For instance: 

V k (x) < 0, Vx6 X . 

V k (x) > \u k l h (x) k \ - 1 Vx fe G X. 

a) General Unreachability and FTT Analysis over Graph Models: The results presented so far in 
this section (Theorem [TJ Corollary [T] and Proposition |4|) are readily applicable to MILMs. These results 



will be applied in Section IV to formulate the verification problem as a convex optimization problem. 
Herein, we present an adaptation of these results to analysis of graph models. 

Definition 5: A cycle C m on a graph G (Af, £) is an ordered list of m triplets (ni, n 2 , kx) , (n 2 , n 3 , k 2 ) , 
{n m , n m+ i, k m ) , where n m+ i = m, and (rij, nj+i, kj) G £, Vj G Z(l,m). A simple cycle is a cycle 
with no strict sub-cycles. 

Corollary 2: Unreachability and FTT Analysis of Graph Models. Consider a program V and its 



graph model G (J\f, 8) . Let V (i, x) = U{ (x) be a Lyapunov invariant for G (Af, S) , satisfying ( 12) and 



<70 (x) < 0, Vx G X (24) 
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and either one of the following two conditions: 

(I) : (Tj (x) > 0, \/x eXiHXi-, i eN\{®} (25) 
(II) : en (x) > 0, VxG X j H T- 1 , z e A/"\ {0} , j G 1 (z) , T G {t*} (26) 

where 

t- 1 = {iG X,\T (x) n ^ 0} 

Then, "P satisfies the unreachability property w.r.t. the collection of sets z G jV\ {0} . In addition, 
if for every simple cycle C G G, we have: 

(9(C)-l)\\a(C)\\ oo + f i(C)>0, and//(C)>0, and ||a (C) |L < oo, (27) 

where 



0(C) = EI ^, MC)= max 4, ||<7 (C)^ = max sup |a,(x)| (28) 

(ij,fc)6C (»j»*)ec (*,.,.)6C xeXAX^ 

then V terminates in at most T u iterations where 

rp sr ((g (c) - 1) Ik (c)!^ + az (c)) - log iiic) v Ha (c)il 

warn + 



Proof: The proof is presented in Appendix II. ■ 
For verification against an overflow violation specified by a diagonal matrix a > 0, Corollary [2] is 



applied with X_ = {1 6 1" ||a > !}• Hence, (25) becomes <Ji(x) > p(x)(\\a 1 x\\ q — 1), 



Va; G Xi, i G Af \ {0} , where p (x) > 0. User-specified assertions, as well as many other standard safety 
specifications, such as absence of division-by-zero can be verified using Corollary [2] (See Table I). 

- Identification of Dead Code: Suppose that we wish to verify that a discrete location z G Af\ {0} in 
a graph model G (Af, S) is unreachable. If a function satisfying the criteria of Corollary [2] with X,^ = M. n 



can be found, then location i can never be reached. Condition (25) then becomes <7j (x) > 0, Vs G M n . 
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TABLE I 

Application of Corollary[2]to the verification of various safety specifications. 





apply Corollary 2 with: 


At location i: 


assert x G X a 


=> 


Xj_ := {x g R" \ x e R n \X a } 


At location i: 


assert x ^ X a 


=> 


XT' f *— TTTi Yt 1 \/' ~\ 

A,_ := {a; € M n | x e A a | 


At location i: 


(expr.)/x 


=> 


Xi_ := {x e R n | x = 0} 


At location i: 


2fc/* r 


=> 


X ( _ := {x G R" | x < 0} 


At location i: 


log (x D ) 


=> 


:= {x eR n \ x D < 0} 


At location i: 


dead code 


=> 


JSQ_ := i?" 



Example 4: Consider the following program 



void ComputcTurnRatc (void) 

L0 : {double x = {0}; double y = {*PtrToY}; 

LI : while (1) 

L2 : { y = *PtrToY; 
L3 : x= (5 * sin(y) + l)/3; 

L4: ifx>-l{ 
L5 : x = x+ 1.0472; 

L6 : TurnRate = y/x; } 

L7 : else { 

L8 : TurnRate = 100 * y/3.1416 }} 




Program 3 



Graph of an abstraction of Program 3 



Note that x can be zero right after the assignment x = (5sin(y) + l)/3. However, at location L6, x 
cannot be zero and division-by-zero will not occur. The graph model of an abstraction of Program 3 
is shown next to the program and is defined by the following elements: T 65 : x i— > x + 1.0472, and 
T41 : x 1 — y [—4/3,2] . The other transition labels are identity. The only non-universal passport labels are 
II54 and n 84 as shown in the figure. Define 

(7 6 (x) = -x 2 - lOOx + 1, 0-5 (x) = -(x + 1309/1250) 2 - lOOx - 2543/25 
<t (x) = G\ (x) = 04 (x) = cr 8 (x) = —x 2 + 2x — 3. 



It can be verified that V (x) = <7j (x) is a (9, l)-Lyapunov invariant for Program 3 with variable rates: 
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6> 65 = 1, and 9^ = for all ^ (6,5). Since 

—2 = sup <t (x) < inf cx 6 (x) = 1 

xGXo xeX - 

the state (6,x = 0) cannot be reached. Hence, a division by zero will never occur. We will show in the 
next section how to find such functions in general. 

IV. Computation of Lyapunov Invariants 
It is well known that the main difficulty in using Lyapunov functions in system analysis is finding 
them. Naturally, using Lyapunov invariants in software analysis inherits the same difficulties. However, 
the recent advances in hardware and software technology, e.g., semi-definite programming [|22|. [|6Ti . 
and linear programming software [|23ll present an opportunity for new approaches to software verification 
based on numerical optimization. 
A. Preliminaries 

1) Convex Parameterization of Lyapunov Invariants: The chances of finding a Lyapunov invariant are 
increased when ^ is only required on a subset of X\Xoo. For instance, for 9 < 1, it is tempting to 
replace Q with 

V (x+) - 9V (x) < -fi, Vx G X\X^ :V(x)<l, x + G / (x) (29) 
In this formulation V is not required to satisfy ([9]) for those states which cannot be reached from X . 



However, the set of all functions V : X i— > K satisfying (29) is not convex and finding a solution for 



<(29j) is typically much harder than ([9]). Such non-convex formulations are not considered in this paper. 



The first step in the search for a function 7:l4R satisfying ([9]) is selecting a finite-dimensional 
linear parameterization of a candidate function V: 

n 

V (x) = V T (x) = nV k (x), r = {r k ) n k=l , r k G R, (30) 

k=l 

where V k '■ X i— > R are fixed basis functions. Next, for every r = {T k ) k= i let 

0(t) = max V T (x + ) — 9V T (x), 
xex\x 00 , x+ef(x) 

(assuming for simplicity that the maximum does exist). Since (•) is a maximum over x of a family of 
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linear functions in r, </>(■) is a convex function. If minimizing </>(■) over the unit disk yields a negative 
minimum, the optimal r* defines a valid Lyapunov invariant V T *(x). Otherwise, no linear combination 
( |30| ) yields a valid solution for 0. 

The success and efficiency of the proposed approach depend on computability of (•) and its subgra- 
dients. While <p (■) is convex, the same does not necessarily hold for V T (x + ) — 6V T {x) as a function of x. 
In fact, if X\Xoo is non-convex, which is often the case even for very simple programs, computation of 
</>(■) becomes a non-convex optimization problem even if V T (x + ) — V T (x) is a nice (e.g. linear or concave 
and smooth) function of x. To get around this hurdle, we propose using convex relaxation techniques 
which essentially lead to computation of a convex upper bound for <p 

2) Convex Relaxation Techniques: Such techniques constitute a broad class of techniques for construct- 
ing finite-dimensional, convex approximations for difficult non-convex optimization problems. Some of 
the results most relevant to the software verification framework presented in this paper can be found 
in 113611 for SDP relaxation of binary integer programs, [l39l and ll47l for SDP relaxation of quadratic 
programs, [|64l for iS-Procedure in robustness analysis, and ll50l . [|49ll for sum-of-squares relaxation in 
polynomial non-negativity verification. We provide a brief overview of the latter two techniques. 

a) The S -Procedure : The 5-Procedure is commonly used for construction of Lyapunov functions 
for nonlinear dynamical systems. Let functions fa : X i->- R, % G Z (0, m) , and ^ : J h> 1, j e Z (1, n) 
be given, and suppose that we are concerned with evaluating the following assertions: 

(I): 0o 0) > 0, Vx e {x e X \ & (x) > 0, if), (x) = 0, « G Z (1, m) , j G Z (1, n)} (31) 

m n 

(II): 3ri G R + , 3(i j G K, such that O (ar) > ^ ^ ( a; ) + Pj&J ^ ' ^ 

i=i j=i 

The implication (II) — > (I) is trivial. The process of replacing assertion (I) by its relaxed version (II) is 
called the iS-Procedure. Note that condition (II) is convex in decision variables Xj and jij. The implication 
(I) — > (II) is generally not true and the 5-Procedure is called lossless for special cases where (I) and (II) 
are equivalent. A well-known such case is when m — 1, n — 0, and 4>q, <pi are quadratic functionals. A 
comprehensive discussion of the 5-Procedure as well as available results on its losslessness can be found 
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in [|26l . Other variations of 5-Procedure with non-strict inequalities exist as well. 

b) Sum-of-Squares (SOS) Relaxation : The SOS relaxation technique can be interpreted as the 
generalized version of the 5-Procedure and is concerned with verification of the following assertion: 

/,■ (x) > 0, Vj E J, g k (x) 7^ 0, Vfc G K, h t (x) =0,VIeL^ -/„ (x) > 0, (33) 

where fj,gk, hi are polynomial functions. It is easy to see that the problem is equivalent to verification 
of emptiness of a semialgebraic set, a necessary and sufficient condition for which is given by the Posi- 
tivstellensatz Theorem In practice, sufficient conditions in the form of nonnegativity of polynomials 
are formulated. The non-negativity conditions are in turn relaxed to SOS conditions. Let S [yi, . . . ,y m ] 

denote the set of SOS polynomials in m variables yi,...,y m , i.e. the set of polynomials that can be 

t 

represented as p = J^Ph Pi ^m, where P m is the polynomial ring of m variables with real coefficients. 



Then, a sufficient condition for (33 ) is that there exist SOS polynomials r , Tj, E £ [x] and polynomials 
pi, such that 

Matlab toolboxes SOSTOOLS (H), or YALMIP J351 automate the process of converting an SOS problem 
to an SDP, which is subsequently solved by available software packages such as LMILAB ll22l . or SeDumi 
[|6T1l . Interested readers are referred to 091, 001, B50l, lEH for more details. 

B. Optimization of Lyapunov Invariants for Mixed-Integer Linear Models 

Natural Lyapunov invariant candidates for MILMs are quadratic and affine functionals. 

I) Quadratic Invariants: The linear parameterization of the space of quadratic functionals mapping 

]R n to K is given by: 



X>1 = { V : W 1 1-> R | V(x) = 



x 
1 



T 

n+1 



P 



P g S" +i } , (34) 



where S n is the set of n-by-n symmetric matrices. We have the following lemma. 
Lemma 1: Consider a program V and its MILM S (F, H, X , n, n w , n v ) . The program admits a quadratic 

(9, p) -Lyapunov invariant V E V%, if there exists a matrix Y E W leXnH , n e = n + n w + n v + 1, a 

diagonal matrix D v E U) nv , a positive semidefinite diagonal matrix D xw E D" +n ™ ; and a symmetric 
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matrix P £ § n+1 , satisfying the following LMIs: 

L\PL X - QL T 2 PL 2 ■< He (YH) + L^D XW L 3 + L{D V L± - (A + //) LjL 5 
A = Trace D xw + Trace D v 

where 















T 




T 




Li = 


F 










, Li = 






0(n e -l)xl 




L 5 




0l x (n e -l) 1 




0(n„+l)x(n+n„) 




. Olxn„ 




1 



Proof: The proof is presented in Appendix II ■ 

The following theorem summarizes our results for verification of absence of overflow and/or FTT for 

MILMs. The result follows from Lemma [T] and Corollary [TJ with q = 2, h = f, though the theorem is 

presented without a detailed proof. 

Theorem 2: Optimization-Based MILM Verification. Let a : -< a < I n be a diagonal positive 

definite matrix specifying the overflow limit. An overflow runtime error does not occur during any 

execution of V if there exist matrices Yj £ l" eX " fl , diagonal matrices D iv E D™ 1 ", positive semidefinite 

diagonal matrices D ixw £ D™ +n ™ ; and a symmetric matrix P £ S> n+1 satisfying the following LMIs: 

[ x Q 1 ]P[ x lf< 0, Vx £ X (35) 
L\PL X - 6L T 2 PL 2 < He (Y 1 H) + L?D lxw L 3 + L T A D lv L A - {\ x + //) L^L 5 (36) 



L'f KL X - L T 2 PL 2 < He (Y 2 H) + L^D 2xw L 3 + L^'D 2l> L4 - A 2 L^ L 5 



where A = diag{a , — 1} , A, = Trace D ixw + Trace An, and ^ A; 



(37) 



1,2. In addition, if 



H + 8 > I, then P terminates in a most T M steps where T u is given in (22). 

2) Affine Invariants: Affine Lyapunov invariants can often establish strong properties, e.g., bound- 
edness, for variables with simple uncoupled dynamics (e.g. counters) at a low computational cost. For 
variables with more complicated dynamics, affine invariants may simply establish sign-invariance (e.g., 
Xi > 0) or more generally, upper or lower bounds on some linear combination of certain variables. As 
we will observe in Section |V} establishing these simple behavioral properties is important as they can 
be recursively added to the model (e.g., the matrix H in a MILM, or the invariant sets Xj in a graph 
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model) to improve the chances of success in proving stronger properties via higher order invariants. The 
linear parameterization of the subspace of linear functionals mapping IR n to R, is given by: 

Vl = \v : R n ^ E | V(x) = K T [x 1] T , K G . (38) 

It is possible to search for the affine invariants via semidefinite programming or linear programming. 

Proposition 5: SDP Characterization of Linear Invariants: There exists a (9, /i)-Lyapunov invariant 
V G for a program V = S (F, H, X ,n,n w ,n v ) , if there exists a matrix F G M" eXn - H , a diagonal 
matrix D v G B n ", a positive semidefinite diagonal matrix D xw G oJ +ni " )x(n+nu,) , and a matrix X G M n+1 
satisfying the following LMI: 

Ee(L^KL 5 - 9LjK T L 2 ) -< Ee(YH) + L%D XW L 3 + LjD v L 4 - (A + /i) L^L 5 (39) 
where A = Trace D XVJ + Trace -Dt, and ^ -D^. 

Proposition 6: LP Characterization of Linear Invariants: There exists a /i)-Lyapunov invariant 
for a program V = S (F, H, X ,n,n wi n v ) in the class V^, if there exists a matrix F G IR lxrtH , and 
nonnegative matrices D v , D v G ]R lxn ", D xw , D xw G ]R lx ( n +M ; and a matrix G M n+1 satisfying: 

K T Li - 6K T L 2 -YH- (D xw - D XW )L 3 - (D v - D V )L, - {D 1 + y) L 5 = (40a) 

D l + (D„ + D v ) l r + (D xw + D xw ) l n+nw < (40b) 

D v , D v , D xw , D xw > (40c) 



where Di is either or —1. As a special case of (40), a subset of all the affine invariants is characterized 



by the set of all solutions of the following system of linear equations: 

K T L X - 6K T L 2 + L 5 = (41) 

Remark 1: When the objective is to establish properties of the form Kx > a for a fixed K, (e.g., when 



establishing sign-invariance for certain variables), matrix K in (39)— (41 ) is fixed and thus one can make 
9 a decision variable subject to 9 > 0. Exploiting this convexity is extremely helpful for successfully 
establishing such properties. 
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The advantage of using SDP for computation of linear invariants is that efficient SDP relaxations for 
treatment of binary variables exists ll39l , ll47l . [|25l . though the computational costs are typically higher 
than the LP-based approaches. In contrast, linear programming relaxations of the binary constraints are 
more involved than the corresponding SDP relaxations. Two extreme remedies can be readily considered. 
The first is to relax the binary constraints and treat the variables as continuous variables. The second is to 
consider each of the 2 Tlv different possibilities (one for each vertex of { — 1, 1}™") separately. This approach 
is practical only if n v is small. More sophisticated schemes can be developed based on hierarchical 
relaxations or convex hull approximations of binary integer programs I1601 , 11361 . 

C. Optimization of Lyapunov Invariants for Graph Models 

A linear parameterization of the subspace of polynomial functional with total degree less than or equal 

to d is given by: 

Vt = iv : W 1 ^ R | V(x) = K T Z (x), K e R N , N = ( U + d \ j (42) 



where Z (x) is a vector of length ( n + ) , consisting of all monomials of degree less than or equal to d in 



n 



variables x±, x n . A linear parametrization of Lyapunov invariants for graph models is defined according 
to (10), where for every i £ Af, we have cr, (■) £ Vt , where d(i) is a selected degree bound for <jj (■) . 



Depending on the dynamics of the model, the degree bounds d {%) , and the convex relaxation technique, 
the corresponding optimization problem will become a linear, semidefinite, or SOS optimization problem. 

1 ) Node-wise Polynomial Invariants: We present generic conditions for verification over graph models 
using SOS programming. For the specific case of verification of linear graph models via quadratic 
invariants and the 5-Procedure for relaxation of non-convex constraints we present explicit LMIs. The 
corresponding Matlab code is available at ll65l . The following theorem follows from Corollary [2} 

Theorem 3: Optimization-Based Graph Model Verification. Consider a program V, and its graph 



model G (Af,£) . Let V : Q n ^ R, be given by (jlOb, where a { (•) £ Vt {i) . Then, the functions a. 



i £ N define a Lyapunov invariant for V, if for all k) £ £ we have: 

- CT 3 -(T* (x,w)) + e%a r (x) - 4 £ E [x,w] subject to (x,w) £ ((X, n nJJ x [-1, n S% (43) 
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Furthermore, V satisfies the unreachability property w.r.t. the collection of sets X { _, i G J\f \ {0} , if there 
exist Ei G (0, oo) , i G J\f\ {0} , such that 



— <T0 (x) G E [x] subject to i 6 
(Tj (x) - Ei G £ [x] subject to x 6 I, fl i G N\ {0} 



(44) 



(45) 



As discussed in Section IV-A2b ( the SOS relaxation techniques can be applied for formulating the search 



problem for functions <Xj satisfying (43)-(45) as a convex optimization problem. For instance, if 

((x, n n* ) x n 5* = {(x, w) \ f p (x, w) > o, h t (x, w) = o], 



then, ( |43| ) can be relaxed as an SOS optimization problem of the following form: 

-°j( T ji (X, W)) + (x) ~ //J -^Tpfp -^Tpgfpfg - J^Mi G £ [x, ttf] , S.t. T p , T p(? G £ [x, w] . 

Software packages such as SOSTOOLS (541 or YALMIP J35) can then be used for formulating the SOS 

optimization problems as semidefinite programs. 

2) Node-wise Quadratic Invariants for Quasi-Linear Graph Models: In a quasi-linear graph model, 

all the transition labels are linear. The passport labels and the invariant sets are defined by linear and/or 

quadratic functions. Let x = [x 1] T , n = n + 1, and let, for all k) G £, a compact description of 



the set Xi n 11^ be given by: 



v, n n* = {x g R n l x T Q k jitl x > o, x T 4 >m x = o, g%x > o, //);,• = 0} 



(46) 



where C ,, G n , ' 

J* 1 



9fi xn ffk rz 



, s G 5*^, and G 



m 



G M|, where S* 



and M fc j are some sets. For all (i,j,k) G £ let the transitions be given by T^x = A^x + B^w + E 



— — — ^ - — - »™ - — - ^ ji 1 j 

where w G [—1, . Before proceeding to the next theorem, for convenience, we define the following 

matrices: 



31 



A 3i 



B 



El 



Olxn 0lxo fc 1 



W 



3* 



I„h O^fc. xl 



Olxn lxo fc 1 



Olxn lx „fc 



(47) 



(48) 
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Theorem 4: SDP-Based Verification of Quasi-Linear Graph Models. Consider a program V, and its 
quasi-linear graph model G (Af, £) . For j G Af, Let ctj y be a diagonal matrix specifying the overflow 
limit and let Aj = diag {cnj 2 , — l} • Let 9^ = for all (0, j, k) G £. An overflow runtime error does not 
occur during any execution of V, if for all (i,j, k) G £ the following LMIs are satisfied: 



T T PjT - B%L T P i L ■< W T DW - (Tr (D) + ^) /C T /C + He (YHC+ZQC) + A£ 



J^kjT 1 ' - C' PiC r< W J 'DW - Tr(D)/C J 7C + He(FUC+Z££) + A 



where for all k) (E £ we have 



ji / ^ ^ ji,s ' / j "ji,m J -' - rL ji,m J ~' 

s m 

Tfc _ \ ^ ~k rT n k n , \ ^ ~k pT -ok n 



(Y, Y, Z, Z, D, D) = (Xi, Y* Z* Z* £>* £*) 



Q, H, K, C, W) = (J?*, G* , H%, K%, L%, W*) 



(49) 
(50) 



(51) 

(52) 

(53) 
(54) 



where the 5-Procedure multipliers in (51)— (53) are constrained as follows: 



V(i,j,k)e£, seS", 



Pji, m , Pji,m e K, V j, fc) g 5, m g Af* 



Yj, y^RWi, V(i,j,k)e£ 

\/{i,j,k)e£ 
V(i,j,k) G £, 



Z fe Z k G M (n+<? ^ +1)x ^ 



(55) 
(56) 
(57) 
(58) 
(59) 



and the matrices in (54) are defined in (46)— (48). In addition, if for every simple cycle C G G, we have: 



e(c) + fi(C)>i, 



v{C) >o, 



(60) 



with (C) and // (C) defined as in (28), then V terminates in at most T u iterations where 



CGG:6»(C)^1 



log (0(C)+/i(C)-l)-log /x(C) 
log 0(C) 



^ /i(C) 



Proof: See Appendix II. 
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Remarks 

1) A simpler, computationally less demanding but possibly more conservative variation of Theorem [4] 



can be formulated by replacing (50) with 

At -Pi^O, Vz G N. (61) 



Subsequently, all the associated multipliers (i.e., the variables denoted by tilde characters in (55 )— (59 )) 
will be eliminated. 



2) It is possible to use the same set of multipliers (i.e., the variables defined in (55)— (59)) across 
all edges (i,j,k) G £. This approach would also result in a computationally less demanding but 
possibly conservative formulation. 

3) To improve the chances of finding stronger invariants, quadratic equality and inequality constraints 
can be generated from linear constraints and added to the set of quadratic constraints in char- 



acterization of the set Xj D (cf. equation 46). Given scalar linear constraints of the form 



G% jr x > 0, r G Z (l,g%) and/or H^ s x = 0, s G Z (l,h%) , define R SlS2 = Hj^ + H^H tl , 



R siri = H 1 Si G ri + G^H Sl , and Q rir2 = G l ri G T2 + G L r2 G rx . Then, we have x 1 R S1 
x T R Siri x = 0, and x T Q Tir2 x > 0, for all r 1: r 2 G Z and si, s 2 G Z (l, /i^) . 



.r 



4) Replacing ( |6T| ) with zAi — Pi -< 0, where z > is a decision variable can improve the scaling 
and numerical conditioning of the associated semidefinite program, resulting in stronger invariants 
and improved analysis. The statement of the theorem on absence of overflow remains unchanged 
under this modification. Only the upper bound on the maximum number of iterations needs to be 
adjusted accordingly since we now have \\a (C)]]^ < z instead of \\a (C)]]^ < 1. The same scaling 



arguments apply to (50). 



5) MATLAB implementations of Theorem [4] and a few variations of it can be found in [|65l . 
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V. Case Studies 

A. Euclidean Division 

In this section we apply the framework to the analysis of Program 4 displayed below. 



/ * EuclideanDivision.c * / 

FO : int IntegerDivision ( int dd, int dr ) 

Fl : {int q = {0}; int r = {dd}; 

F2 : while (r >= dr) { 
F3: q = q+l; 

F4 : r = r — dr; 

Fix : return r; } 

L0 : int main ( int X, int Y ) { 

LI : int rem = {0}; 

L2 : while (Y > 0) { 

L3 : rem = IntegerDivision (X , Y); 

L4 : X = Y; 

L5 : Y = rem; 

Lm : return X; }} 



Program 4: Euclidean Division and its Graph Model 

Program 4 takes two positive integers X £ [1,M] and Ye [1,M] as the input and returns their greatest 

common divisor by implementing the Euclidean Division algorithm. Note that the Main function in 

Program 4 uses the IntegerDivision program (Program 1). 

1 ) Global Analysis: A global model can be constructed by embedding the dynamics of the Inte- 
gerDivision program within the dynamics of Main. A labeled graph model is shown alongside the 

text of the program. This model has a state space X —J\f x [— M, M] , where M is the set of nodes as 

shown in the graph, and the global state x = [X, Y, rem, dd, dr, q, r] is an element of the hypercube 

[— M, M] . A reduced graph model can be obtained by combining the effects of consecutive transitions 

and relabeling the reduced graph model accordingly. While analysis of the full graph model is possible, 

working with a reduced model is computationally advantageous. Furthermore, mapping the properties of 

the reduced graph model to the original model is algorithmic. Interested readers may consult [|59l for 

further elaboration on this topic. For the graph model of Program 4, a reduced model can be obtained 
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by first eliminating nodes F M , L 4 , L 5 , L 3 , F , Fi, F 3 , F 4 , and Li, (Figure [T] Left) and composing the 
transition and passport labels. Node L 2 can be eliminated as well to obtain a further reduced model with 
only three nodes: F 2 , L , L^. (Figure [T] Right). This is the model that we will analyze. The passport and 



-@ 




^"F2F2 




n F 2 F 2 r F 2 F 2 



Fig. 1. Two reduced models of the graph model of Program 4. 

transition labels associated with the reduced model are as follows: 

— 1 — 2 

T F2F2 : xk[X, Y, rem, dd, dr, q + 1, r - dr] T F2F2 : x^[Y, r, r, Y, r, 0, Y] 



T L0F2 :x4[X, Y, 0, X, Y, 0, X] 



n F2F2 : {x | 1 < r < dr - 1} II F2F2 : {x | r > dr} 



7f2Ln : x^[Y, r, r, dd, dr, q, r] 
n F2 L* : {x | r < dr - 1, r < 0} 



Finally, the invariant sets that can be readily included in the graph model (cf. Section II-B2a) are: 

X L0 = {x | M > X, M > Y, X > 1, Y > 1} , X F2 = {x | dd = X, dr = Y} , X Lk = {x | Y < 0} . 

We are interested in generating certificates of termination and absence of overflow. First, by recursively 
searching for linear invariants we are able to establish simple lower bounds on all variables in just two 
rounds (the properties established in the each round are added to the model and the next round of search 
begins). For instance, the property X > 1 is established only after Y > 1 is established. These results, 
which were obtained by applying the first part of Theorem [3] (equations (43)-(44) only) with linear 
functional are summarized in Table II. 

TABLE II 



Property 


q > 


Y > 1 


dr > 1 


rem > 


dd > 1 


X > 1 


r > 


Proven in Round 


I 


I 


I 


I 


II 


II 


II 


°"F 2 (x) = 


-q 


1 - Y 


1 - dr 


—rem 


1 - dd 


1 - x 


— r 


(^F2F2' / i F2F2j 


(1,1) 


(1,0) 


(1,0) 


(1,0) 


(0,0) 


(0,0) 


(0,0) 


(^F2F2' ^F2F2) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 
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We then add these properties to the node invariant sets to obtain stronger invariants that certify FTT and 
boundedness of all variables in [— M, M]. By applying Theorem [3] and SOS programming using YALMIP 
|35]|, the following invariants are founcj^after post-processing, rounding the coefficients, and reverifying): 



a 1F2 (x) = 0.4 (Y - M) (2 + M - r) 

\2 



a 2F2 (x) = (q x Y + r) 2 - M 2 
a 3F2 (x) = (q + r) z - M 2 ct 4 f 2 (x) = 0.1 (Y - M + 5Y x M + Y 2 - 6M 2 ) 

a 5F2 (x) = Y + r - 2M + Y x M - M 2 a 6Fa (x) = r x Y + Y - M 2 - M 

The properties proven by these invariants are summarized in the Table III. The specifications that the 
program terminates and that x G [— M, M] 7 for all initial conditions X, Yg [1, M] , could not be established 
in one shot, at least when trying polynomials of degree d < 4. For instance, cr 1F2 certifies boundedness 
of all the variables except q, while a 2 F 2 an d o"3F 2 which certify boundedness of all variables including q 
do not certify FTT. Furthermore, boundedness of some of the variables is established in round II, relying 
on boundedness properties proven in round I. Given a (x) < (which is found in round I), second 
round verification can be done by searching for a strictly positive polynomial p (x) and a nonnegative 
polynomial q (x) > satisfying: 

q (x) a (x) - p (x) ((TxX 2 -M 2 )>0, Te{Tl 2F2 



where the inequality ( [62] ) is further subject to boundedness properties established in round I, as well as 
the usual passport conditions and basic invariant set conditions. 



TABLE III 



Invariant ctf 2 ( x ) = 


CTlF 2 ( x ) 


0"2F 2 ( x ) > °"3F 2 (x) 


°"4F 2 (x) 


05F 2 (x) , cr 6 p 2 (x) 


(^F2F2' ^F2F2) 


(1,0) 


(1,0) 


(1,0) 


(1,1) 


(^F2F2' ^F2F2) 


(1,0.8) 


(0,0) 


(1,0.7) 


(1,1) 


Round I: x 2 < M 2 for x i = 


Y, X, r, dr, rem, dd 


q, Y, dr, rem 


Y, X, r, dr, rem, dd 


Y, dr, rem 


Round II: x 2 < M 2 for x; = 




X,r,dd 




X,r,dd 


Certificate for FTT 


NO 


NO 


NO 


YES, T u = 2M 2 



2 Different choices of polynomial degrees for the Lyapunov invariant function and the multipliers, as well as different choices for 8, /i 
and different rounding schemes lead to different invariants. Note that rounding is not essential. 



EXTENDED VERSION, AUGUST 2011 31 

In conclusion, cr 2 F 2 ( x ) or cr 3F 2 ( x ) m conjunction with cr 5F2 (x) or cr 6F2 (x) prove finite-time termination 
of the algorithm, as well as boundedness of all variables within [— M, M] for all initial conditions X, Y e 
[1, M] , for any M > 1. The provable bound on the number of iterations certified by ct 5 f 2 (x) and cr 6 p 2 (x) 
is T u = 2M 2 (Corollary [2]). If we settle for more conservative specifications, e.g., x G [— kM,kM] 7 for 
all initial conditions X, Y £ [1,M] and sufficiently large k, then it is possible to prove the properties in 
one shot. We show this in the next section. 

2) MIL-GH Model: For comparison, we also constructed the MIL-GH model associated with the 
reduced graph in Figure [TJ The corresponding matrices are omitted for brevity, but details of the model 
along with executable Matlab verification codes can be found in [65 J. The verification theorem used 
in this analysis is an extension of Theorem [2] to analysis of MIL-GHM for specific numerical values 
of M, though it is certainly possible to perform this modeling and analysis exercise for parametric 
bounded values of M. The analysis using the MIL-GHM is in general more conservative than SOS 
optimization over the graph model presented earlier. This can be attributed to the type of relaxations 
proposed (similar to those used in Lemma [T]) for analysis of MILMs and MIL-GHMs. The benefits 
are simplified analysis at a typically much less computational cost. The certificate obtained in this way 
is a single quadratic function (for each numerical value of M), establishing a bound 7 (M) satisfying 
7 (M) > (X 2 + Y 2 + rem 2 + dd 2 + dr 2 + q 2 + r 2 ) 1 ^ 2 . Table IV summarizes the results of this analysis 
which were performed using both Sedumi 13 and LMILAB solvers. 



TABLE IV 



M 


10 2 


10 3 


10 4 


10 5 


10 6 


Solver: LMILAB (22): 7 (M) 


5.99M 


6.34M 


6.43M 


6.49M 


7.05M 


Solver: SEDUMI (61]: 7 (M) 


6.00M 


6.34M 


6.44M 


6.49M 


NAN 


(^F2F2> ^F2F2) 


(1,10- 3 ) 


(1,10-3) 


(1,10-3) 


(1,10-3) 


(I.IO- 3 ) 


(^F2F2> ^¥2¥2) 


(1,10- 3 ) 


(1,10-3) 


(1,10-3) 


(1,10-3) 


(1,10-3) 


Upperbound on iterations 


T u = 2e4 


T u = 8e4 


T u = 8e5 


T u = 1.5e7 


T u = 3e9 



3) Modular Analysis: The preceding results were obtained by analysis of a global model which was 
constructed by embedding the internal dynamics of the program's functions within the global dynamics of 
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the Main function. In contrast, the idea in modular analysis is to model software as the interconnection of 
the program's "building blocks" or "modules", i.e., functions that interact via a set of global variables. The 
dynamics of the functions are then abstracted via Input/Output behavioral models, typically constituting 
equality and/or inequality constraints relating the input and output variables. In our framework, the 
invariant sets of the terminal nodes of the modules (e.g., the set X K associated with the terminal node Y K 
in Program 4) provide such I/O model. Thus, richer characterization of the invariant sets of the terminal 
nodes of the modules are desirable. Correctness of each sub-module must be established separately, while 
correctness of the entire program will be established by verifying the unreachability and termination 
properties w.r.t. the global variables, as well as verifying that a terminal global state will be reached in 
finite-time. This way, the program variables that are private to each function are abstracted away from 
the global dynamics. This approach has the potential to greatly simplify the analysis and improve the 
scalability of the proposed framework as analysis of large size computer programs is undertaken. In this 
section, we apply the framework to modular analysis of Program 4. Detailed analysis of the advantages 
in terms of improving scalability, and the limitations in terms of conservatism the analysis is an important 
and interesting direction of future research. 

The first step is to establish correctness of the IntegerDivision module, for which we obtain 

a 7F2 (dd, dr, q, r) = (q + r) 2 - M 2 

The function a 7F2 is a (1, 0) -invariant proving boundedness of the state variables of IntegerDivision. 
Subject to boundedness, we obtain the function 

<7 8F2 (dd, dr, q, r) = 2r — llq - 6Z 

which is a (1, l)-invariant proving termination of INTEGERDIVISION. The invariant set of node F M can 
thus be characterized by 

Xy, = {(dd, dr, q, r) G [0,M] 4 |r<dr-l} 

The next step is construction of a global model. Given X M , the assignment at L3: 

L3 : rem = IntegerDivision (X , Y) 
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can be abstracted by 

rem = W, s.t. W G [0, M] , W < Y — 1, 
allowing for construction of a global model with variables X, Y, and rem, and an external state-dependent 
input W G [0, M] , satisfying W < Y — 1. Finally, the last step is analysis of the global model. We obtain 
the function <7gL2 (X, Y, rem) = Y x M — M 2 , which is (1, l)-invariant proving both FTT and boundedness 
of all variables within [M, M] . 

B. Filtering 

Program 4 has been adopted with some minor modifications from a similar program analyzed by 
ASTREE JS), Il66ll . The only modification that we have made is the addition of the real-time input 
w G [—1,1] in line L4. Note that when B = (in line L4) we have the same program reported in ll66l . 
In Program 4, the function saturate(.) truncates the real-time input *PtrToInput so that w G [— 1, 1] is 
guaranteed. We first build a graph model of this program. The variables of the graph model are: 

Z, Y, E[0], E[l], S[0], S[l], INIT. 

The variable INIT is Boolean which we model by v G { — 1,1} in the following way: 

INIT = True v = 1, and INIT = False <=> v = -1. 
The graph model of this program is shown in Figure [2j 
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FO 

Fl 

f2 
f3 
f4 
F5 
f6 
F7 
F8 
f9 
FlO 

Fll 
F12 
F13 

LO 

Ll 

L2 
L3 

L4 

L5 

L6 
L7 



/* filter.c */ 

typedef enum {FALSE = 0, TRUE = 1} BOOLEAN; 

BOOLEAN INIT; float Y={0}, Z={0}; 

void filter () { 

static float E[2], S[2]; 

if (INIT) { 

S[0] = Z; 

Y = Z; 
E[0] = Z; 

} else { 

Y = (((((0.5*Z) - (E[0]*0.7)) + (E[l]*0.4)) + (S[0]*1.5)) - (S[l]*0.7)); 
} 

E[l] = E[0]; 
E[0] = Z; 
S[l] = S[0]; 
S[0] = Y; 

} 

void main () { 
Z = 0.2 * Z + 5; 
INIT = TRUE; 
while (1) { 

wait (0.001); w = saturate(*PtrToInput); /*updates real-time input*/ 
Z = 0.9 * Z + 35+B*w; 
filter (); 

INIT = FALSE; 
} 

} 



Program 4: Example of filtering in safety-critical software 
Example adapted from reference ll66ll with minor modifications. 



In order to construct a more compact model, several lines of code corresponding to consecutive 
transitions have been combined, and only the first line corresponding to a series of consecutive transitions 
has been assigned a node on the graph. Specifically, the combined lines are: {Ll, L2}, and {F3, F4, F5} , 
and {F9, FlO, Fll, F12} . The net effect of the eliminated lines is captured in the model by taking the 
composition of the functions that are coded at these lines. 
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Fig. 2. The graph of Program 4. 

Letting x= [Z, Y, E[0], E[l], S[0], S[l], v] T , the state space of the graph model is TV X IR 7 , where 
TV = {L0,Ll,L3,L4,L5,L6,Lx} U {FO, F2, F3, F6, F7, F9, F13} , as shown in Figure g The initial 
condition is fixed: X® := {LO} x x where x = [0, 0, 0, 0, 0, 0, 1] T . The corresponding passport and 
transition labels are as follows: 



21.3L1 


: x^ [0.2Z + 5, Y, E[0], E[l], S[0], S[l], 1] 








: x^ [0.9Z + 35 + Bw, Y, E[0], E[l], S[0], 


S[l], v] 




21.3L6 


: x^[Z, Y, E[0], E[l], S[0], S[l], -1] 






^F9F3 


: x^[Z, Z, Z, E[l], Z, S[l], v] 






Tf9F7 


: x i— > [Z, 0.5Z - 0.7E[0] + 0.4E[1] + 1.5S[0] - 


0.7S[1], E[0], E[l], 


S[0], S[l], v 


2"f13F9 


:x^[Z, Y, Z, E[0], Y, S[0], v] 






riF3F2 


= {x | v = 1} , n F6F2 = {x | v = -1} , 


riL4L3 = K i 


HlmL3 = 



The other transitions are identity maps. Theorem[4]can be applied for verification of absence of overflow 
(note that the program does not terminate). We apply Theorem [4] for two different values of B, B = 
which corresponds to no real-time input, and B = 20. For each value of B G {0, 20}, absence of overflow 
is verified using two different methods. In Method I, we set A 3 ■ = A = diag{Mr 2 l 7 , —1} , and thus, we 



EXTENDED VERSION, AUGUST 2011 



36 



obtain a certificate for M > ||x|| 2 using one quadratic invariant. In Method II, for each i G {1, .., 7} we set 
Aj = A = diag {M~ 2 ej, — 1} , where is the i-th. standard unit vector in IR 7 (cf. Corollary [I] and Equation 



<|23j)). As a result, we obtain a much tighter upperbound in the form of M > at the expense of more 

computations. Table V summarizes the results of this analysis. Analysis of the floating-point operations 
is based on the abstractions discussed in Appendix I. The corresponding Matlab codes can be found 
in [|65l . In this analysis, the solver SeDuMi version 1.3 lloTI is used for semidefinite optimization, and 
Yalmip 113511 for compilation. For comparison we point out that the static analyzer Astree reports a bound 
of M = 1327.05 for the B = case 11661 . 

TABLE V 





Computations: 


infinite-precision 


64-bit double precision 


32-bit single precision 




Method I 


B = 


M = 884.95 


M = 884.96 


M = 884.96 


M > ||x|| 2 > HxlU 


Method II 


B = 


M = 373.11 


M = 373.12 


M = 373.12 




Method I 


B = 20 


M = 1400.95 


M = 1402.81 


M = 1402.87 


M> ||x|| 2 > HxlU 


Method II 


B = 20 


M = 609.83 


M = 609.83 


M = 609.83 


M > Hex, 




(fij.Mij) 


(0.98,0) 


(0.98,0) 


(0.98,0) 





VI. Concluding Remarks 
We took a systems -theoretic approach to software analysis, and presented a framework based on convex 
optimization of Lyapunov invariants for verification of a range of important specifications for software 
systems, including finite-time termination and absence of run-time errors such as overflow, out-of-bounds 
array indexing, division-by-zero, and user-defined program assertions. The verification problem is reduced 
to solving a numerical optimization problem, which when feasible, results in a certificate for the desired 
specification. The novelty of the framework, and consequently, the main contributions of this paper are in 
the systematic transfer of Lyapunov functions and the associated computational techniques from control 
systems to software analysis. The presented work can be extended in several directions. These include 
understanding the limitations of modular analysis of programs, perturbation analysis of the Lyapunov 
certificates to quantify robustness with respect to round-off errors, extension to systems with software in 
closed loop with hardware, and adaptation of the framework to specific classes of software. 
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VII. Appendix I 

A. Semialgebraic Set-Valued Abstractions of Commonly-Used Nonlinearities 
- Trigonometric Functions: 

Abstraction of trigonometric functions can be obtained by approximation of the Taylor series expan- 
sion followed by representation of the absolute error by a static bounded uncertainty. For instance, an 
abstraction of the sin (•) function can be constructed as follows: 



Abstraction of sin (x) 


x <= \—7L ZLl 

L 2 ' 2 J 


X £ [— TT, 71"] 


sin (x) € {x + aw | w G [—1, 1]} 


a = 0.571 


a = 3.142 


sin (x) £ {x — ^x 3 + aw \ w £ [—1, 1]} 


a = 0.076 


a = 2.027 



Abstraction of cos (•) is similar. It is also possible to obtain piecewise linear abstractions by first 
approximating the function by a piece-wise linear (PWL) function and then representing the absolute 
error by a bounded uncertainty. Section II-B (Proposition [T]) establishes universality of representation 
of generic PWL functions via binary and continuous variables and an algorithmic construction is given 
in the proof of the proposition in Appendix II. For instance, if x G [0, 7r/2] then a piecewise linear 
approximation with absolute error less than 0.06 can be constructed in the following way: 

S={(x,v,w) |a; = 0.2[(l + ^)(l + ^ 2 ) + (l-^)(3 + ^ 2 )], (w,v) G [-1, l] 2 x {-1,1}} (63a) 
sm 0) e {Tx E | x E e S} , T : x E \-> 0.45 (l + v)x + (l-v) (0.2a; + 0.2) + 0.06w 1 (63b) 



- The Sign Function (sgn) and the Absolute Value Function (abs): 

The sign function (sgn(x) = llro oo) ( x ) ~ H(-oo,o) ( x )) ma y appear explicitly as one of the program's 
functions, or implicitly as a model for conditional switching. A particular abstraction of sgn(-) is as 
follows: sgn(x) G {v \ xv > 0, v E { — 1, 1}}. Note that sgn (0) is equal to 1, while the abstraction is 
multi-valued at zero: sgn(0) G { — 1, 1} . The absolute value function can be represented (precisely) over 
[—1, 1] in the following way: 

abs (x) = {xv | x = 0.5 (v + w) , (w, v) G [—1, 1] x { — 1, 1}} 

- Modulo Arithmetic: 
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Consider the function mod :ZxZ->Z defined in the following way: 

mod (t, s) = t — ns, where n = \_-\ . 
Abstraction of the function mod(.,.) over bounded subsets of the set of integers is as follows. Let 

t G [0, Ms), where M > is an integer. Then: 

mod(M) G \ t-^sJ2( 1 + v k) I (t-ks)v k > 0, v k e {-1,1}, k = 1, M - 1 1 . (64) 



k=l 



It follows from (64) that the result of summation modulo s of two integers ti,t 2 G [0, Ms) can be 



abstracted in the following way: 

2M-1 

h + t 2 - -s ( X + I (*i + *2 - M^fc > 0, w fe G {-1, 1} , fc = 1, 2M - 1 
fc=i 

B. Floating-Point and Fixed-Point Arithmetic 

This subsection is based on fl42), Chapter 7. Interested readers are referred to fl42), [|43l for a more 
comprehensive discussion of computations with floats. For computations with floating-point numbers, the 
IEEE 754-1985 norm has become the hardware standard, and is supported by most popular programming 

languages such as C. In this standard, a float number is represented by a triplet (s, /, e) , where: 

• s G {0, 1} is the sign bit. 

• / is the fractional part, represented by a p-bit unsigned integer: f\. . . f p , fi G {0, 1}. 

• e is the biased exponent, represented by a g-bit unsigned integer: e\ . . . e q , e« G {0, 1}. 
A floating point number z = (s, f, e) is then in one of the following forms: 

. z = (_l) s x T~ Uas x 1./, when 1 < e < 2 q - 2. 

. z = (-l) s x 2 l - bias x 0./, when e = 0, and / ^ 0. 

• z = +0 (when s — 1) and z = — (when s — 0) and e = / = 0. 

• z = +oo (when s = 1) and z = — oo (when s = 0) and e = 2 9 — 1, / ^ 0. 
. z = iVaiV when e = 2 9 - 1, / = 0. 

The values of p, q, and bias depend on the specific format: 

• If format is 32-bit single precision (f=32), then bias = 127, q — 8, and p = 23. 

• If format is 64-bit double precision (f=64), then bias = 1023, q = 11, and p = 52. 
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Other formatting standards such as long double or quadruple precision also exist. In floating point 
computations, the result of performing the arithmetic operation © G {+,— , x, -=-} on two float variables 
x and y is stored in afloat variable z := float (x © y, f) which is a complicated function of x, y, and the 
format f. Examples of the format f include the IEEE 754 with 32-bit single precision, or IEEE 754 with 
64-bit double precision. A floating-point operation is equivalent to performing the operation over the set 
of real numbers followed by rounding the result to afloat ll43l. Il42l. In IEEE 754 the possible rounding 
modes are rounding towards 0, towards +oo, towards — oo, and to the nearest (n). The rounding function 
r f m : M. — > FU {e} maps a real number to a float number or to runtime error e, depending on the format 
'f and the rounding mode 'm'. We refer the reader to [42J for more details on the rounding function. 
For our purposes, it is sufficient to say that for all rounding modes, the following relation holds: 

\/x G [— CK f , CK f ] : |r f>m (x) — x I < 7f \x\ + (5 f 

where j f := 2 _p , and atf := (2 — 2~ p ) 2 29 ~ bms ~ 2 is the largest non-infinite number, and /? f := 2 1 ~ bms ~ p is 

the smallest non-zero positive number. 

Based on the above discussions, an abstraction of the floating-point arithmetic operators can be 
constructed in the following way: 



x + y g 


[—CKf, CKf] = 


» float 


(x + y) 


= z G {x 


+ y + 5w 


W G 


[-1,1], ^ 


= Tf (\ x \ 


+ M) + M 


x — y G 


[-CKf, CKf] = 


» float 


(x-y) 


= z G {x 


— y + Sw 


w G 


[-1,1], ^ 


= 7f (kl 


+ H)+/5f} 


x x y G 


[-CKf, CKf] = 


» float 


(x X y) 


= z G {x 


x y + 5w 


w G 


[-1,1], 5 


= 7f (N 


x M) + A-} 


x 4- y G 


[-CKf, CKf] = 


> float 


(x + y) 


= z G {x 


-i- y + 6w 


W G 


[-1,1], 5 


= 7f (N 


-M) + /9f} 



where the constants «f , /3f and 7f are defined as before. The above abstractions can still be complicated 
as the magnitude of 5 depends on the operands x and y. If all of the program variables (including the 
result of the arithmetic operation) reside in [—a, a] where a CKf then a simpler but more conservative 
abstraction can be constructed in the following way: 

x © y G [—a, a] => float (x ® y) — z G {x ® y + 5w \ w G [— 1, 1] , 5 = a^f + /3f} 
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For instance, if f=32 and a = 10 6 , then 5 = 0.12, and if f=64 and a = 10 10 , then 5 = 2.3 x 1(T 6 . 
Abstractions of arithmetic operations in fixed-point computations is similar to the above. The magnitude 
of 5 will depend on the number of bits and the dynamic range. For instance, in the two's complement 
format we have: 5 = p (2 b — l) 1 , where p is the dynamic range and b is the number of bits. A case 
study in software verification based on the above discussed abstractions of floating-point computations 



has been presented in Section V-B 



C. Graph models with state-dependent or time-varying edges 

Consider the following fragment from a program with two variables: i6l, and i G {1, 2} 



LI : x = A [i] x; 
L2 : expression 



Then, the state transition from node 1 to 2 is defined by T 2 \ : (x, i) — > (A [ % ] x, i) , where A is an array 
of size 2. An algorithm similar to the one used in construction of MILMs (cf. Proposition^ can be used 
to construct a transition label using semi-algebraic set-valued maps. For the above example, this can be 
done in the following way: 

T 2 i(x,i) = {T 21 (x,i,v 1 ,v 2 ) | (ar,i,«i,v 2 ) G S 21 }. 
T 21 (x,i,vt,v 2 ) = A [1] v\x + A [2] v 2 x. 

S 21 = {(x,i,v 1 ,v 2 ) | v x + v 2 = 1, v u v 2 G {0, 1} , v x (i - 1) + v 2 (i - 2) = 0} . 

In light of Proposition [TJ generalization of this technique to multidimensional arrays with arbitrary finite 
size is straightforward. However, the number of binary decision variables grows (linearly) with the number 
of elements in the array, which can be undesirable for large arrays. An alternative approach to constructing 
graph models with fixed labels is to attempt to derive a fixed map by computing the net effect of several 
lines of code IBTI . When applicable, the result is a fixed map which is obtained by taking the composition 
of several functions. This is an instance of a more general concept in computer science, namely, extracting 
the higher level semantics, which allows one to define more compact models of the software. 
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&K>1 




k + 1 

(x,y) ^ T tl [k](x,y) 



{k<N}\ 



{k = N} 



(x,y) i-> r. 2 [fc](x,y) 



(*,>-) k> r.,p,/|(«,j/) 





7^1 









{i = N,j = N} 



(x,y) ^ T, 5 [i,j](x,y) 



i l-» 1+1 

0,y) ^ r. 5 p,7](x,^) 



Fig. 3. Graph model of a code fragment (Program 5) with state-dependent labels. The transition labels are shown in boxes and the passport 
labels are in brackets. For simplicity, only the non-identity transition labels are shown. 



For instance, consider the following code fragment: 



LI : 


for ( k = 1 ; k == N ; k++) { 




L2 : 


y[ k ] = x[ k ]; x[ k ] = 0; } 




L3 : 


for ( i = 1 ; i == N ; i++) { 




L4 : 


for ( j = 1 ; j == N ; j++) 


{ 


L5 : 


x[i] =x[i ] +y[j ] * 


A[i][j ]; 


L6 : 


}} 





Program 5: A code fragment from a linear filtering application 



A graph model of this code fragment is shown in Figure |3j where: 

y=[y[l]---y[N)) T , x = [x [1] ■ • • x [N)f 



T* 2 [k] : 



T* 5 ■ 



y 












-> 
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X 



EXTENDED VERSION, AUGUST 2011 



-12 



where is an N x iV matrix which is zero everywhere except at the (i,j)-th entry which is 1, and A t j 
denotes the (z, j)-th entry of A (Aij = A [i] [j]) which is an V x A" array in the code lETTl . The remaining 
transition labels correspond to the counter variables i, j, k, and have been specified on the graph in 
Figure [3j From node 1 to node 3, the net effect is: 



N 

UT.2 [k] 
k=l 



/ 




and from node 3 to node 6 the net effect is: 



N N 

nn^M 

i=l,j=l 



/ 

A 



VIII. Appendix II 



Proof of Proposition 



N 



The proof is by construction. Let X = [jXi, where the ATj's are compact 

i=l 

polytopic sets. Further, assume that each X{ is characterized by a finite set of linear inequality constraints: 

Xi = {x \ < Si, Si G R N * Xn , Si eR N '}. Let v G {-1, 1}^ , and consider the following sets: 



N 



(x,v) 



J2>Vi = -N + 2, (Six - Si) (vi + 1) < 0, Vi G {-1,1} : % G Z(1,JV) 



i=l 
N 



G xy = \{x,y) | ^(1 + Vi) (AiX + Bi) = y, (x, v) G G 2 
i=i 



First, we prove that G xy = G\. Define N binary vectors: rf G { — 1,1}^, j G Z(1,AT), according to 
the following rule: ^ = 1 if and only if i = j. Also define I (x) = {j G Z (1, N) \ x G Xj} . Now, let 

(x ,/(x )) G Gi. Then 

(x ,^) G G xv : Vj G I(a;o) . 

Therefore, (x 2AjX + 25j) G G^j, : Vj G I (x ) , which by supposition, implies that (x , / (x )) G G xy . 
This proves that G\ C G^. Now, let (xo,yo) G G^. Then, there exists v G {—1,1}, such that 
(x ,v) G G m . Hence, there exists j G Z(l, N) , such that x G X J; and y = 2Ajx + 2Bj. It follows 
from the definition of / that (xo,2/o) £ Gi. This proves that G xy C Gi. We have shown that Gi = G xy . 
It remains to show that G xy has a representation equivalent to G2. 
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Define new variables m = xv,i E [— 1, l} n . Then 



X,V) I V = Yh=1 ( A i X + A i U i + B i + BiVi) , SiUi + SiX - SiVi - Si < 



fi A 
^xy — 



(65) 



Eti Vi = -N + 2, m = xv i} Vi E {-1, 1} : i E Z (1, iV) 



The nonlinear map (x,t>i) H- ttj can be represented by an affine transformation involving auxiliary 
variables z { E [—1, l] n , and Zj G [—1, l] n , subject to a set of linear constraints, in the following way: 



(66) 



equivalently: 



Ui = 2Zi - X - (Vi - 1) 1„, Zi = X- Zi- l n , Z { = (V{ - 1) 1„ + W i5 = (-Uj - 1) l n + 117*, 

where w^Wi E [— l,l] n . Since by assumption each Xi is bounded, for all i E Z(l,iV) and all j E 
Z (1, Ni) , the quantities 



R iH = min Sj,x — S;,- 



(67) 



exist, are finite and can be computed by solving iVjX N linear programs given in (67 ) (Sj,- and s« denote 



the j-th row of Sj and s» respectively). Let R { = diag {i? i7 }. Then G xy as defined in (65 ) is equivalent 



to: 



G 



X ,V) \ V = J2f=l ( A i X + A i U i + B i + BiVi) , (X, Ui, Vi) E H 



(68) 



where H.—{(x, V{) \ x,Ui,Vi satisfy (69) for some z^z^w^Wi E [— 1, l] n , and u>i E [—1,1]^} 



H= < 



= S^i + Six - s^i - Si- Ri (wi + Ijv. 
Ui = 2z { - x - (vi - 1) l n 

Z { = (Vi - 1) l n + 



-X Zj l n , 2j, Zj 6 [ 1, 1] 
■■(-Vi - 1) l n + Wi, 



> (69) 



Eii^i =-iV + 2, ^€{-1,1} 



The equality constraints that define the Matrix H are precisely the equalities in (69), while the equality 



constraint in (68) defines the Matrix F. This completes the proof. 



Proof of Proposition [2| First, consider the unreachability property. Assume to the contrary, that V 
does not satisfy unreachability w.r.t. X— Then, there exists a solution X = x(.) of S(X, f, Xq, Xoo), 
and a positive integer t- E Z + such that x (0) E X , and x G X_. It follows from X C X , and 
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f(x) f(x) that X = x{.) is also a solution of S(X, f ^Xq^X^). Therefore, we have: 

x 6L,LCl x (*_) G X_. 

However, x(t) G X_ contradicts the fact that S(X, f, X , Xqo) satisfies safety w.r.t. X-. Proof of 
FTT is similar: let X = x(.) be any solution of <S(X, /, X^X,*,). Since X = x(.) is also a solution of 
<S(X, f, X , Xoo), it follows that there exists t T G Z + such that x (t T ) G I M . Since x(t T ) is also an 
element of X, it follows that x (tx) G Im fl X. Since by definition PI X C holds, we must have 
z (*r) e ^oo- ■ 



Proof of Proposition^ Note that ( 13 )— (fl5|) imply that V is negative-definite along the trajectories 



of iS, except possibly for V (x (0)) which can be zero when 77 = 0. Let X be any solution of S. Since 
V is uniformly bounded on X, we have: 

-mi oo <y(x(i))<0, Vx(i)G*, t>i. 

Now, assume that there exists a sequence X = (x(0), x(l), . . . , x(t), . . . ) of elements from X satisfying 
([I]), but not reaching a terminal state in finite time. That is, x (t) Vi G Z+. Then, it can be 



verified that if t > T u , where T u is given by (16), we must have: V (x (£)) < — ||V|| , which contradicts 
boundedness of V. ■ 

Proof of Theorem [!]■ Assume that 5 has a solution X = (x (0) , ...,x (£_) , ...) , where x (0) G X and 
x (t_) G X_. Let 

7h = inf V (x) 

cce/i- 1 ^-) 

First, we claim that j h < max{V r (x (£_)), V(x (t_ — 1))} . If /i = /, we have x (t_) G (X_) and 
7h < V(x (t-)). If h — f, we have x (t_ - 1) G /tT 1 (X_) and 7^ < V(x (t_ - 1)), hence the claim. 
Now, consider the 9 = 1 case. Since V is monotonically decreasing along solutions of S, we must have: 

7 h = inf V(x) < max{V(x(t_)), V(x(t_ - 1))} < V(x(0)) < supK(x) (70) 



which contradicts (17). Note that if \x > and h = I, then (70) holds as a strict inequality and we can 
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replace (17) with its non-strict version. Next, consider case (I) , for which, V need not be monotonic 



along the trajectories. Partition X into two subsets X and X_ Q such that X = X U X and 

V (x) < Vx G X , and V (x) > ViG X 

Now, assume that S has a solution (x (0) x (t_ ),...) , where x (0) G X and x (£_) G X_. Since 
V (x (0)) > and 9 < 1, we have 1/ (x (t)) < V (x (0)) , Vt > 0. Therefore, 

j h = inf V(x) < max{V(x(t-)), V(x(t- - 1))} < V(x(0)) < supK(x) 



which contradicts (17). Next, assume that S has a solution Xj= (x (0) , ...,x(t_) , ...) , where x (0) G X 



and x (t_) G X_. In this case, regardless of the value of 0, we must have V (x (i)) < 0, Vt, implying that 



7/i < 0, and hence, contradicting (fT8J). Note that if /i = I and either /i > 0, or 9 > 0, then (18) can be 



replaced with its non-strict version. Finally, consider case (II). Due to (19), V is strictly monotonically 



decreasing along the solutions of S. The rest of the argument is similar to the 9 = 1 case. 



Proof of Corollary U\ It follows from (21 ) and the definition of X_ that: 



V(x) > sup|||a" 1 /i(x)|| g - l} > sup{||a" 1 /i(x)|| oo - 1} > 0, Vx G X. (71) 



It then follows from (71 ) and (20) that: 



inf V (x) > > sup V(x) 

Hence, the first statement of the Corollary follows from Theorem [T] The upperbound on the number of 
iterations follows from Proposition |4] and the fact that swp xeX ^ x _ uXoo y \V (x)| < 1. ■ 

Proof of Corollary |2j- The unreachability property follows directly from Theorem [T] The finite time 
termination property holds because it follows from (12), (24) and (27) along with Proposition |4} that the 



maximum number of iterations around every simple cycle C is finite. The upperbound on the number of 
iterations is the sum of the maximum number of iterations over every simple cycle. ■ 

Proof of Lemma [7j Define x e = (x,w,v, 1) T , where x G [—1, l] n , w G [—1, l] nw , v G { — 1, l}™ 1 " . 
Recall that (x, 1) T = L 2 x e , and that for all x e satisfying Hx e = 0, there holds: (x + , 1) = (Fx e , 1) = L±x e . 
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It follows from Proposition [3] that Q holds if: 

x T e L T x PL x x e - 9x T e LlPL 2 x e < -p,, s.t. Hx e = 0, L 3 x e E [-1, l] n+Uw , L A x e E {-1, l} n " . (72) 



Recall from the iS-Procedure ((31 ) and (32)) that the assertion a (y) < 0, E [— 1, l] n holds if there exist 
nonnegative constants n > 0, i = 1, n, such that a (y) < (yf — 1) = y T ry — Trace (r) , where 
r = diagjrj} ^ 0. Similarly, the assertion a (y) < 0,V?/ E { — 1,1}™ holds if there exist constants pi 
such that cr (y) < Yl Pi (Hi ~ 1) = V T PV ~ Trace (p) , where p = diag {p^} . Applying these relaxations 



to (72), we obtain sufficient conditions for (72) to hold 



L T x PL x x e -Qx T e L\PL 2 x e < x T e (YH + H T Y T ) x e +x^LlD xw L 3 x e +x^LlD v L 4 x e -p-Tmce{D xw +D v ) 



Together with ^ D xw , the above condition is equivalent to the LMIs in Lemma [T] ■ 
Proof of Theorem^ Since = for all (0, j, k) E £, we have <jj (x) < 0, for all j in the outgoing 
set of node 0. This replaces condition (24) in Corollary [2j The first part of the theorem thus follows from 
Corollary [2] and an application of the iS-Procedure relaxation technique in a similar fashion to the proof 
of Lemma [T] The second part of the theorem also follows from Corollary [2] by noting that due to (50), 
we have ||er (C)|| < 1 for every simple cycle C EG. ■ 
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